#include <assert.h>
#include <stdio.h>
#include <stdlib.h>

#include "../src/em.c"

typedef struct _EmTestDatum EmTestDatum;
struct _EmTestDatum 
{ 
  struct
    {
      float     kappa;
      float     logalpha;
      float     tau[(5-1) * 4];
      float     loglambda[4];
      float     priorz[5];
      int       clamp;
      Rating    ratings[8];
      float     logalphastddev;
    }                   inputs;

  struct
    {
      float    newq;
      float    logalpha;
      float    pz[5];
    }                   desired_output;
};

static void
test_em (void)
{
  EmTestDatum data[] = 
    {{{0.252243,  0.874555, {0.418849, 0.9538, 0.962466, 0.864147, 0.860376,  0.112554, 0.0449485, 0.514772, 0.213757, 0.477305, 0.907866,  0.366047, 0.423551, 0.379591, 0.55147, 0.0902069}, {0.143739,  0.839714, 0.328195, 0.738698}, {0.382411, 0.942393, 0.0759522,  0.864143, 0.963562},  3, {{1, 1}, {3, 1}, {2, 3}, {0, 4}, {3, 2}, {3, 2}, {3, 2}, {0,  2}}, 1.11349}, {-14.007080991750735170, -0.92989546710448087080,  {6.16678e-11, 8.58471e-8, 4.35379e-6, 0.999996,  1.80479e-8}}}, {{0.999995,  0.103186, {0.876039, 0.0685383, 0.485224, 0.889429, 0.398734,  0.160672, 0.119176, 0.465878, 0.0191431, 0.609202, 0.0289696,  0.32214, 0.179429, 0.281007, 0.290272, 0.939728}, {0.237037,  0.205055, 0.426129, 0.976166}, {0.248443, 0.0915681, 0.426134,  0.87298,  0.372404}, -1, {{0, 0}, {0, 3}, {2, 4}, {2, 0}, {0, 4}, {0,  0}, {2, 0}, {1, 0}},  1.94091}, {-15.283687501422783161, -2.4553607389522542022,  {0.290298, 0.0865795, 0.267953, 0.298494, 0.0566753}}}, {{0.983551,  0.973671, {0.862358, 0.821734, 0.517672, 0.954528, 0.253156,  0.792764, 0.195533, 0.775098, 0.972148, 0.502493, 0.255805,  0.538062, 0.767094, 0.0763634, 0.279638, 0.289618}, {0.675525,  0.65023, 0.406658, 0.917214}, {0.652496, 0.709319, 0.423108,  0.943543, 0.790138},  1, {{1, 3}, {1, 2}, {2, 4}, {0, 1}, {3, 3}, {1, 2}, {3, 2}, {2,  3}}, 1.90544}, {-11.328312367516183977, -0.22777377875266204369,  {6.84175e-22, 0.0000750046, 0.00834876, 0.990545,  0.00103102}}}, {{0.989016,  0.536982, {0.0948208, 0.709903, 0.213918, 0.564834, 0.592328,  0.454098, 0.675856, 0.79774, 0.515965, 0.17446, 0.386238,  0.122215, 0.865735, 0.767801, 0.469024, 0.469719}, {0.156416,  0.344693, 0.52548, 0.679582}, {0.268831, 0.439258, 0.536465,  0.142599, 0.17401},  3, {{3, 3}, {2, 3}, {1, 2}, {3, 2}, {0, 0}, {3, 1}, {0, 0}, {2,  0}}, 1.32255}, {-16.914192078423761047, -1.5373231002829225221,  {1.95227e-6, 8.07763e-6, 9.58222e-6, 0.99998,  1.42195e-7}}}, {{0.577766,  0.581682, {0.275257, 0.646691, 0.780025, 0.0657171, 0.100797,  0.260453, 0.65781, 0.199982, 0.332996, 0.79143, 0.188091,  0.0435655, 0.988303, 0.265949, 0.50851, 0.774735}, {0.549045,  0.729485, 0.36591, 0.600724}, {0.81969, 0.406937, 0.788145,  0.0190424, 0.544433},  4, {{2, 2}, {3, 0}, {1, 4}, {2, 3}, {3, 1}, {0, 1}, {3, 0}, {2,  1}}, 1.00812}, {-19.343135926948070635, -2.1132031077260230622,  {0.000022317, 0.0000298721, 0.0000233613, 0.0000102134,  0.999914}}}, {{0.953325,  0.443636, {0.499793, 0.350309, 0.753344, 0.11064, 0.708363,  0.162218, 0.709778, 0.122337, 0.442414, 0.653708, 0.935044,  0.573292, 0.71293, 0.287798, 0.334319, 0.753602}, {0.305992,  0.499653, 0.315277, 0.20917}, {0.545746, 0.491534, 0.361951,  0.765534,  0.0459527}, -1, {{2, 0}, {3, 1}, {0, 2}, {3, 1}, {1, 1}, {1,  3}, {0, 2}, {3, 2}},  1.60861}, {-12.310134053909681944, -0.28135987054007140711,  {9.96366e-6, 0.0510565, 0.9179, 0.0310331,  2.10203e-8}}}, {{0.654894,  0.337589, {0.979008, 0.89883, 0.532557, 0.895175, 0.3253, 0.963786,  0.959265, 0.182245, 0.037502, 0.629467, 0.205662, 0.876253,  0.537849, 0.31419, 0.996493, 0.330507}, {0.0463144, 0.952239,  0.230959, 0.284555}, {0.905089, 0.343631, 0.576065, 0.946965,  0.926081}, -1, {{0, 4}, {3, 3}, {1, 2}, {1, 0}, {2, 1}, {2,  3}, {0, 3}, {3, 1}},  1.04351}, {-15.100407085477124001, -1.1331153103297776505,  {0.000954905, 0.0233148, 0.396163, 0.551764,  0.0278036}}}, {{0.0517903,  0.600781, {0.481015, 0.0842432, 0.869545, 0.563279, 0.851548,  0.878581, 0.993292, 0.0254307, 0.537358, 0.882088, 0.662784,  0.979116, 0.585119, 0.651129, 0.37823, 0.0740273}, {0.241488,  0.075064, 0.431264, 0.147946}, {0.796686, 0.0315561, 0.379474,  0.547165, 0.315671},  1, {{2, 1}, {1, 1}, {2, 2}, {2, 0}, {0, 1}, {0, 3}, {1, 1}, {0,  1}}, 1.50993}, {-12.180690029775833389, -0.92171997699058773059,  {3.50022e-7, 0.999997, 2.37443e-6, 1.58426e-7,  6.44233e-10}}}} ;

  for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
    {
      ParameterizedNormalDistribution prior_logalpha = 
        parameterized_normal_distribution (data[i].inputs.logalphastddev);

      float logalpha = data[i].inputs.logalpha;
      float pz[5];
      float logpriorz[5];

      for (unsigned int j = 0; j < 5; ++j)
        {
          pz[j] = drand48 ();
          logpriorz[j] = log (data[i].inputs.priorz[j]);
        }

      float q = em (data[i].inputs.kappa,
                    &logalpha,
                    5,
                    data[i].inputs.tau,
                    data[i].inputs.loglambda,
                    pz,
                    logpriorz,
                    data[i].inputs.clamp,
                    data[i].inputs.ratings,
                    8,
                    &prior_logalpha.base,
                    1000,
                    1e-4);

      assert (fabs (q - data[i].desired_output.newq) <=
              1e-3 * (1 + fabs (q) + fabs (data[i].desired_output.newq)) ||
              ( fprintf (stderr, "%u: %g ?= %g\n", i, q, data[i].desired_output.newq),
                0 ));
      
      assert (fabs (logalpha - data[i].desired_output.logalpha) <=
              1e-3 * (1 + fabs (logalpha) + fabs (data[i].desired_output.logalpha)) ||
              ( fprintf (stderr, "(oops logalpha ... but this is ok) %g ?= %g\n", logalpha, data[i].desired_output.logalpha),
                1 ));

      for (unsigned int j = 0; j < 5; ++j)
        {
          assert (fabs (pz[j] - data[i].desired_output.pz[j]) <=
                  1e-3 * (1 + fabs (pz[j]) + fabs (data[i].desired_output.pz[j])) ||
                  ( fprintf (stderr, "(oops pz[%u] ... but this is ok) %g ?= %g\n", j, pz[j], data[i].desired_output.pz[j]),
                    1 ));
        }
    }
}

int 
main (void)
{
  srand48 (69);

  test_em ();

  return 0;
}
